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Abstract. The two planets about the star GJ 876 appear to have undergone extensive migration from their point 
of origin in the protoplanetary disk — both because of their close proximity to the star (30 and 60 day orbital 
periods) and because of their occupying three stable orbital resonances at the 2:1 mean-motion commensurability. 
The resonances were most likely established by converging differential migration of the planets leading to capture 
into the resonances. A problem with this scenario is that continued migration of the system while it is trapped in 
the resonances leads to orbital eccentricities that rapidly exceed the observational upper limits of ei « 0.31 and 
62 ~ 0.05. As seen in forced 3-body simulations, these lower eccentricities would persist during migration only for 
an eccentricity damping rate 62/62 exceeding « 40 < 12 / 02 . Previous theoretical and numerical analyses have found 
e/e ~ <i/a or even eccentricity growth through disk-planet interactions. 

In an attempt to find effects that could relax the excessive eccentricity damping requirement, we explore the 
evolution of the GJ 876 system using two-dimensional hydrodynamical simulations that include viscous heating 
and radiative cooling in some cases. Before we evolve the whole system, the disk with just the outer planet 
embedded is brought into equilibrium. We find that the relaxed disk remains circular in all models for low 
pfanet-to-star mass ratios 52 , but becomes eccentric for high mass ratios for those models with fixed temperature 
structure. The disk in models with full radiative thermodynamics remains circular for all </2 considered due to the 
larger disk temperatures. Given the small stellar mass, the mass ratio for the GJ 876 system is rather high (with 
minimum <72 = 5.65 x 10“®), and so the GJ 876 disk may have been slightly eccentric during the migration. 
With a range of parameter values, we find that a hydrodynamic evolution within the resonance, where only 
the outer planet interacts with the disk, always rapidly leads to large values of eccentricities that exceed those 
observed — similar to the three-body results. The resonance corresponding to the resonant angle 81 = 2 A 2 — 
Ai — tui (involving the inner planet’s periapse longitude, tui) is always captured first. There is no additional 
delay in capturing 82 = 2 A 2 — Ai — -^[72 into resonance that is attributable to the secular prograde contribution 
to the precession of zi 72 from the interaction with the disk, but an eccentric disk can induce a large outer planet 
eccentricity 62 before capture and thereby further delay capture of 82 for larger planetary masses. The delay 
in capturing 82 into libration, while delaying the resonance-induced growth of 62 , has no effect on the forced 
eccentricities of both planets, which are uniquely determined by the resonance conditions, once both 8 j are 
librating. 

Only if mass is removed from the disk on a time scale of the order of the migration time scale (before there 
has been extensive migration after capture), as might occur for photoevaporation in the late phases of planet 
formation, can we end up with eccentricities that are consistent with the observations. 
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1. Introduction 

Among the 14 known extrasolar planetary systems with 
multiple planets, at least three exhibit pairs of planets that 
are likely to be in orbital resonances at low order commen- 
surabilities of their mean motions. The pair of planets 

Send offprint requests to: W. Kley, 
e-mail: kleyOtat .physik.uni-tuebingen.de 


around GJ 876 (Marcy et al. 2001) is well confirmed to 
be deep in the resonances associated with the 2:1 mean 
motion commensurability (Laughlin & Chambers 2001; 
Rivera & Lissauer 2001; Laughlin et al. 2004), where res¬ 
onant angles 0i = 2 X 2 — Ai — vui, 02 = 2 X 2 — Ai — - 072 , 
and Aot = W 2 — wi = 0i — 02 are all librating about 0° 
with small amplitudes. Here Xj are mean longitudes and 
vjj are longitudes of periapse, both numbered from the 
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inside out. The pair of planets around HD 82943 is likely 
to also be in resonances at 2:1 (Mayor et al. 2004), and 
the middle pair of planets in the 4 planet system orbiting 
55 Cnc may be in resonances at 3:1 (Marcy et al. 2002; 
McArthur et al. 2004). Thus as many as one-fourth of 
multiple-planet systems contain planets in mean motion 
resonances. 

Both theoretical and numerical analyses verify the 
ubiquity of planet migration due to interaction with the 
protoplanetary disk of gas and dust (e.g. Ward 1997; 
Nelson et al. 2000). The clearing of disk material between 
two planets both capable of opening a gap in the disk 
leads to differential migration of the two planets as at 
least the outer planet is forced in by the material re¬ 
maining outside its orbit. The convergence of the orbits 
naturally allows capture into stable resonant orbital con¬ 
figurations. Hydrodynamical simulations of two embedded 
planets capturing each other into resonance have been per¬ 
formed by several groups (Kley 2000; Bryden et al. 2000; 
Snellgrove et al. 2001; Papaloizou 2003; Kley et al. 2004). 
The resonance capture has also been analyzed with exten¬ 
sive three-body calculations, where migration is simply 
imposed with either the semi-major axis migration rate 
a/a and the eccentricity damping rate e/e specified ex¬ 
plicitly or ad hoc forces added to produce migration and 
eccentricity damping (Snellgrove et al. 2001; Lee & Peale 
2002; Nelson & Papaloizou 2002; Kley et al. 2004; Lee 
2004). 

Depending on the migration rates, masses, and initial 
orbital separations and eccentricities of the two planets, 
capture can occur in different resonances. For planets as 
massive as those in GJ 876, capture into the 2:1 resonances 
is robust if the initial a^jax ^ 2 and the initial eccentrici¬ 
ties are small, and it is most likely that the 2:1 resonances 
for GJ 876 were established by converging differential mi¬ 
gration (Lee & Peale 2002). The sequence of 2:1 resonance 
configurations that a system with initially nearly circu¬ 
lar orbits is driven through by continued migration de¬ 
pends mainly on the planetary mass ratio Mi/M 2 . For 
M 1 /M 2 « 0.31 as in the GJ 876 system, a system is first 
captured into antisymmetric configurations with 61 librat- 
ing about 0 ° and Atu (and hence 62 ) librating about 180°. 
Gontinued migration forces ei to larger values and 62 from 
increasing to decreasing until 62 ~ 0 when ei « 0.1 Then 
the system converts to symmetric configurations like that 
of GJ 876, with both 0i and Aw librating about 0° (Lee 
2004; see also Fig. 16 below). There are other types of 2:1 
resonance configurations (with both 61 and Aw librating 
about 180° or asymmetric librations of 0i and Aw about 
values far from either 0° or 180°) for M 1 /M 2 « 0.31, but 
they are either unstable for planets as massive as those 
in GJ 876 or not reachable by convergent migration of 
planets with nearly constant masses and coplanar orbits 
(Lee 2004). Asymmetric libration configurations can re¬ 
sult from convergent migration for 2:1 resonances with 
M 1 /M 2 ^ 0.95 and for 3:1 resonances for a wider range 
of M 1 /M 2 (Beauge et al. 2003; Ferraz-Mello et al. 2003; 
Kley et al. 2004; Lee 2004). 


While it is easy to understand the symmetric resonance 
configuration of GJ 876 from convergent migration, the 
small eccentricities of this system are a puzzle. The con¬ 
tinued migration after capture into resonance drives the 
pair of planets deeper into the resonances ( 2 n 2 — ni < 0 
increasing, with Uj = Xj being the mean orbital motions). 
The identical mean retrograde precessions of the Wj must 
therefore decrease in magnitude to keep 0 j near zero, i.e., 
to maintain the resonance configuration. This decrease is 
effected by increasing the orbital eccentricities when the 
system is in the configuration with both 9i and 62 librat¬ 
ing about 0°. The three-body calculations have shown that 
the eccentricity growth is rapid if there is no eccentricity 
damping, with the average eccentricities ((ei) = 0.255, 
( 62 ) = 0.035) of the sinf = 0.78 fit of Laughlin & 
Chambers (2001) exceeded after only a 7% decrease in 
the orbital radii after capture (Lee & Peale 2002). The ec¬ 
centricities can be maintained at the observed low values 
as migration continues if there is significant eccentricity 
damping of e/e « lOOa/a (Lee & Peale 2002). In contrast 
the full hydrodynamical calculations typically give simi¬ 
lar timescales for both migration and eccentricity damp¬ 
ing (Kley et al. 2004). The low eccentricities become even 
more problematic if disk-planet interactions drive eccen¬ 
tricity growth. Ogilvie & Lubow (2003) and Goldreich & 
Sari (2003) have suggested that the eccentricity-damping 
corotation torques can be reduced sufficiently to trigger 
eccentricity growth if the eccentricity is above a critical 
value, and the eccentricity of the outer planet in GJ 876 
is well above their estimates for this critical value. 

Hydrodynamical simulations of two planets interacting 
with their protoplanetary disk are computationally expen¬ 
sive, and the number of such simulations in the literature is 
relatively small (Kley 2000; Bryden et al. 2000; Snellgrove 
et al. 2001; Papaloizou 2003; Kley et al. 2004). In addition, 
often the masses used in these simulations are not appro¬ 
priate for the GJ 876 system or the simulations did not 
continue for a very long time. While there has been more 
extensive three-body calculations with imposed migration 
and eccentricity damping (Snellgrove et al. 2001; Lee & 
Peale 2002; Nelson & Papaloizou 2002; Kley et al. 2004; 
Lee 2004), they do not model disk-planet interactions self- 
consistently and, in particular, the effects of the apsidal 
precession induced by the disk have not been considered. 

In this paper we present a more comprehensive set of 
numerical calculations treating the evolution of two plan¬ 
ets interacting with their protoplanetary disk, with spe¬ 
cial emphasis on the resonant system GJ 876. For this 
purpose we solve the full hydrodynamical equations, in¬ 
cluding viscous heating and radiative cooling effects in 
some cases, and follow the joint planet-disk evolution over 
several thousand orbits of the planets. We find that we 
cannot reproduce the observed low eccentricities using 
these straightforward assumptions. Taking into account 
mass loss from the disk, as might be the case in the late 
dissipation stages of protoplanetary disks subject to, e.g., 



Kley, et al.: Modeling GJ 876 


3 


Table 1. Orbital parameters of the two planets of the plan¬ 
etary system GJ 876 at epoch JD 2449679.6316 assuming co¬ 
planarity and i = 90° as given by Laughlin et al. (2004). The 
adopted stellar mass is M* = O.32M0. P denotes the orbital 
period, M the mass of the planet, a the semi-major axis, e the 
eccentricity, w the angle of periastron at epoch, and q the mass 
ratio M/Mt. 



P 

M 

a 

e 


Q 


[d] 

[Mjup] 

[AU] 


[deg] 

[10-^] 

Inner 

30.38 

0.597 

0.13 

0.218 

154 

1.78 

Outer 

60.93 

1.89 

0.21 

0.029 

149 

5.65 


photoevaporation, we are able to reduce the final values 
of the eccentricities close to the observed values. 

In Sect. 2 we give an overview of the observational data 
for GJ 876. In Sect. 3 we explain our physical and numer¬ 
ical approach, and in Sect. 4 we present the results of our 
simulations. This is followed in Sect. 5 with three-body 
calculations and analytic theory to interpret the results of 
Sect. 4. We state our conclusions in Sect. 6. 

2. The System GJ 876 

The first planet orbiting GJ 876 was discovered in 1998 
(Marcy et al. 1998; Delfosse et al. 1998). With the discov¬ 
ery of the second inner planet (Marcy et al. 2001), the near 
2:1 commensurability of the orbital periods (« 30 and 60 
days) implied stable mean-motion resonances, which were 
soon confirmed. The combined minimum mass of the plan¬ 
ets is about 0.0074 of the stellar mass (M* = O. 32 M 0 ). 
The large relative planetary masses and short orbital pe¬ 
riods meant that a dynamical fit to the radial velocity 
data that accounted for the mutual gravitational interac¬ 
tion of the planets was required. This was accomplished 
by Laughlin & Ghambers (2001) and Rivera & Lissauer 
(2001), who used a Levenberg-Marquardt minimization 
scheme driving an iV-body integrator to find best-fit ini¬ 
tial orbital elements. Recently new data obtained with the 
Keck Telescope have also been included in the analysis of 
the system (Laughlin et al. 2004), and the latest best-fit 
orbital parameters are shown in Table 1 for an assumed 
coplanar orbital inclination of f = 90°. 

Although a dynamical fit can in principle yield the or¬ 
bital inclinations and masses of the planets, the value 
of the fit for GJ 876 is relatively insensitive to the incli¬ 
nation i of the assumed coplanar orbits for 35° ^ i < 90° 
and rises rapidly only for i ^ 35° (Laughlin et al. 2004). In 
Fig. 1 we show the planetary mass ratio M 1 /M 2 and the 
average orbital eccentricities (ei) and ( 62 ) for the best-fits 
found by Laughlin et al. (2004) for i > 35°. The deviation 
from Keplerian motion that is most strongly constrained 
by the available observations of GJ 876 is the resonance 
induced retrograde apsidal precession of the orbits at an 
average rate oiw = —41° yr“^, and indeed this precession 
has now been observed for more than one full period. As 
i decreases, the individual planetary mass Mj increases 
roughly as 1 /sini (and Mi/M 2 is nearly constant and 
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Fig. 1. Planetary mass ratio. Mi/M 2 , average orbital eccen¬ 
tricities, (ei) and (62), and the ratio of eccentricity damping 
rate to migration rate, K, for the GJ 876 best-fit solutions of 
Laughlin et al. (2004) with coplanar inclination i > 35°. The 
value of A = |e2/e2|/[a2/a2| is for the equilibrium eccentrici¬ 
ties to match (cj) if there is inward migration and damping of 
the outer planet only. 


« 0.31). The increase in planetary mass acts to increase 
w. This increase can be offset by an increase in (e^ ), which 
acts to decrease zu, keeping rb near the observed value 
without any significant change in until i ^ 35°. To al¬ 
low for the mass uncertainties, we shall consider systems 
with different planetary masses in our numerical models. 

The most striking feature of GJ 876 is the exact 2:1 
orbital resonance of the planets. The angles di, 62 , and 
Aw describe the state of the system for coplanar orbits, 
and all three angles librate about 0° in the GJ 876 system. 
The libration amplitudes are |0i|max = 7°, |02|max = 34°, 
and |A-ci7|i„ax = 34° for the i = 90° fit shown in Table 1, 
and the smallest libration amplitudes in the range 35° < 
i < 90° are |0i|max « 5° at i « 50° and |02|max « 16° 
and |Ati7|i„ax ~ 13° at i « 35° (Laughlin et al. 2004). The 
small libration amplitudes indicate that the system is deep 
in the resonance, which occurs naturally in the differential 
migration scheme for forming the resonant structure if the 
planets approach the resonance slowly with initially small 
eccentricities (Lee & Peale 2002). 

Lee & Peale (2002) have shown using three-body inte¬ 
grations with imposed migration and eccentricity damp¬ 
ing that the eccentricities reach equilibrium values that 
remain constant for arbitrarily long migration within the 
resonances if e^/ej = —K\dj/aj\. For inward migration 




4 


Kley, et al.: Modeling GJ 876 


and damping of the outer planet only, K « 100 is re¬ 
quired to produce equilibrium eccentricities that match 
the eccentricities of the sini = 0.78 (or i = 51°) fit of 
Laughlin & Chambers (2001). We have repeated the cal¬ 
culations in Lee & Peale (2002) for the updated fits in 
Laughlin et al. (2004), and the value of Lf as a function of 
i is shown in Fig. 1. For i = 51°, K « 80, which is close to 
the value inferred by Lee & Peale (2002). Benedict et al. 
(2002) have claimed from HST astrometric measurements 
that the inclination is very close to 90°, which would re¬ 
quire a significantly larger amount of eccentricity damping 
{K « 170). The increase in {ej) with decreasing i leads to 
a decrease in K, but K must still be greater than about 
40 if i ^ 35°. 

3. The Hydrodynamical Model 

Even though the physical setup is different from exist¬ 
ing simulations, the models presented here are calculated 
roughly in the same manner as those described previously 
in Kley (1998, 1999) for single planets and in Kley (2000) 
for multiple planets. The reader is referred to those papers 
for details on the computational aspects of this type of 
simulations. Other similar models, following explicitly the 
motion of single and multiple planets in disks, have been 
presented by Nelson et al. (2000), Bryden et al. (2000), 
Snellgrove et al. (2001), Papaloizou (2003), and more re¬ 
cently for resonant configurations by Kley et al. (2004). 

We use cylindrical coordinates (r, ip, z) and consider 
a vertically averaged, infinitesimally thin disk located at 
z = 0, where the origin of the coordinate system is at the 
position of the star. The basic hydrodynamic equations 
(mass and momentum conservation) describing the time 
evolution of such a two-dimensional disk with embedded 
planets have been stated frequently and are not repeated 
here (see Kley 1999). Additional information on the treat¬ 
ment of embedded planets is given in Kley et al. (2004). 

3.1. Energy Equation 

The majority of models presented here use a fixed temper¬ 
ature distribution which follows from the assumption of a 
constant ratio of vertical height H (r) to radial distance r 
from the star. Here we assume H{r)/r = const. = 0.05, 
from which T{r) oc r~^ follows. In this case the pressure 
p is given hy p = Ecg, where S is the surface density 
and Cs the isothermal sound speed. In the situation of a 
given iL/r-ratio there is no need for solving an extra en¬ 
ergy equation, and we refer to those models as isothermal 
(even though the radial temperature is varying). 

We also present radiative models with an improved 
thermodynamic treatment using the thermal energy equa¬ 
tion 

dUc T 

+V ■ (ScvTu) = -pV u + D-Q (1) 

Here, u = (ur,u^p) is the two-dimensional velocity, T the 
(midplane) temperature of the disk, Cv the ratio of specific 


heats, D the dissipation function, and Q is the local ra¬ 
diative cooling. The form of the dissipation function D in 
cylindrical coordinates can be found in Mihalas & Weibel 
Mihalas (1984), and to calculate the radiative losses (from 
the two sides of the disk) we follow D’Angelo et al. (2003a) 
and Gunther et al. (2004) and write 

Q = ( 2 ) 

where ctr is the Stefan-Boltzmann constant and Teff is an 
estimate for the effective temperature (Hubeny 1990) 

Teff Teff = with Teff = ^T+^ + ^ (3) 

For a two-dimensional disk we approximate the mean ver¬ 
tical optical depth by 

T = (4) 

where for the Rosseland mean opacity k the analytical 
formulae by Lin & Papaloizou (1985) are adopted. 

In addition to the simplified treatment of radiation as 
expressed in (1) we also have considered models including 
radiative diffusion in the (r, (p)-plane. In this case a term 

XrdnT^ 

-2HV-F with F = -VT (5) 

pn 

has been added to the rhs of Eq. (1). Here F denotes the 
radiative flux in the (r, (p)-plane, c is the speed of light, a 
the radiation constant, p = S/(2i7) the midplane density, 
and A the flux-limiter (see Kley 1989). 

We work in a rotating reference system, rotating ap¬ 
proximately with the initial period of the outer planet. As 
the coordinate system is accelerated and rotating, we take 
care to include the indirect terms. 

3.2. Initial Setup 

The two-dimensional (r — tp) computational domain con¬ 
sists of a complete ring of the protoplanetary disk cen¬ 
tered on the star. By previous simulations of two embed¬ 
ded planets interacting with the protoplanetary disk it has 
been shown that the inner part of the disk (inside of the 
owter planet) clears rapidly. The final configuration is such 
that two planets orbit the star inside a cavity (Kley et al. 
2004). In this state only the outer planet is still in touch 
with the disk, see Fig. 1 in Kley et al. (2004). Taking ac¬ 
count of this fact, we simplify the setup and choose the 
radial extent of the computational domain (ranging from 
Tmin to Tmax) such that the inner planet orbits entirely 
inside r^nin- This reduces the necessary radial range cov¬ 
ered and leads to a significant reduction in computational 
effort. A typical example of the grid structure is displayed 
below in Fig. 2. In the azimuthal direction for a complete 
annulus we have pmin = 0 to (/3max = 27r. 

The initial hydrodynamic structure of the disk (den¬ 
sity, temperature, velocity) is axisymmetric with respect 
to the location of the star, and the surface density scales 
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as S(r) = with a superimposed initial gap 

(Kley 2000). The initial velocity is pure Keplerian rotation 
(ur — 0,u^ = (GM^/ry/^), and the initial temperature 
stratification is given by T(r) (x r~^ which follows from 
an assumed constant H/r. For the isothermal models the 
initial temperature profile is left unchanged, while for the 
radiative cases it is evolved according to Eq. (1). 

The kinematic viscosity is parameterized by an a- 
description v = aCgH, where the isothermal sound-speed 
is given by Cg = (cvT)^^^, and H{r) is either held fixed 
(for the isothermal models) or, for the radiative models, 
calculated from the temperature/sound-speed as H{r) = 
Cs / nK(c), where 



is the Keplerian angular velocity of the disk. 

3.3. Boundary conditions 

In an effort to ensure a uniform environment for all mod¬ 
els and minimize disturbances (wave reflections) from the 
outer boundary, we impose at rmax damping boundary 
conditions where the density and both velocity compo¬ 
nents are relaxed towards their initial state as 



X [AU] 


Fig. 2. The global grid structure with gray-shaded initial sur¬ 
face density superimposed. Every second grid points is shown. 
The dots denote the location of the star and the two planets 
and the oval line refers to the Roche lobe of the outer planet. 


dX 

dt 


X-X{t = 0) 
rilamp 


Rir? 


( 6 ) 


where X G {T,,Ur,u^}, ruamp = l/f^K(?'max) and i?(r) 
is a dimensionless linear ramp-function rising from 0 to 
1 within [rdampj ?'max]- This damping setup is defined in 
more detail in an international comparison test project^. 
As the initial radial velocity is vanishing, this damping 
routine ensures that no mass flows through the outer 
boundary at rmax- However, in some models described be¬ 
low, matter is allowed to leave the outer boundary. 

At the inner radial boundary rmin outflow conditions 
are applied; matter may flow out, but none is allowed to 
enter. This procedure mimics the accretion process onto 
the star. The density gradient is set to zero at rnCm, while 
the angular velocity there is fixed to be Keplerian. In the 
azimuthal direction, periodic boundary conditions for all 
variables are imposed. 

These specified boundary conditions allow for a well 
defined quasi-stationary state if the planets are not al¬ 
lowed to respond to the disk. 

3.4. Model parameters 

The computational grid in all models has the same radial 
extent from rmin = 0.25 to rmax = 1.20AU. It is covered by 
111 X 450 {Nr X N^) gridcells which are spaced logarith¬ 
mically in radius and equidistant in azimuth. The radius 
beyond which the damping procedure defined above grad¬ 
ually sets in is given by rdamp = IT- The stellar mass 


is O.32M0 in all models. For the planetary masses differ¬ 
ent choices have been made as the exact parameters for 
GJ 876 are not known. The majority of models assumed 
an edge on system i = 90° which leads to (see Tab. 1) 
mi = 0.597Mjup and m 2 = 1.89Mjnp. Recall that the 
index 1 refers to the inner and the index 2 to the outer 
planet. In terms of their mass ratios (M/M*) these values 
are equivalent to qi = 1.78 x 10“^ and q 2 = 5.65 x 10“^. 
In the analysis of the models we prefer to state the plane¬ 
tary masses in terms of mass ratios rather than the actual 
values. Most of the mass ratios used in this investigation 
are close to the edge-on case. The two planets are placed 
initially at semi-major axes of ai = 0.20 AU and 02 = 0.35 
AU. Both planets have an initial eccentricity of 0.01 which 
is comparable to what is found typically in numerical sim¬ 
ulations of disk-planet interactions. Even if the planets 
were started with zero initial eccentricity at those aj, the 
dynamical interaction between the planets alone would 
also lead to non-zero eccentricities of the same order. 

After an initial relaxation procedure, where the orbital 
parameters of the planets are not changed and the disk 
structure is brought to an equilibrium, the planets are ‘re¬ 
leased’ in all cases and are allowed to migrate (change their 
semi-major axes) in accordance with the gravitational disk 
forces/torques exerted on them. 

During the evolution, material may enter from the disk 
into the Roche lobe of a planet. This material is partly 
removed from the simulation to model accretion onto the 
planet; for the detailed procedure see Kley (1999). In some 


http://www.astro.su.se/'pawel/planets/test.hydro.htmlnodels this mass is not added to the dynamical mass of 
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the planets, so that mass is not strictly conserved. In other 
simulations the mass is added to that of the planet, main¬ 
taining conservation of mass. In Table 2 below it is in¬ 
dicated that only one of the presented models (h8a) has 
a varying dynamical mass. In general it does not change 
the outcome of a simulation noticably whether mass is 
accreted onto a planet or not. 

For the viscosity a value of a = 0.01 is used for all mod¬ 
els. This is probably on the large side for protoplanetary 
disks, but it allows for a rapid evolution of the system and 
hence a reasonable computational effort; a larger a speeds 
up the evolution and migration of the planet. It has been 
shown earlier that the migration speed has no influence on 
final magnitude of the eccentricities (Lee & Peale 2002) 
which tends to justify this approach. In addition, a larger 
a-value leads to a gap that is not so well cleared (Kley 
1999) which will tend to induce a larger periastron ad¬ 
vance (d?) and an increased eccentricity damping. Both 
effects presumably serve to minimize the final eccentricity 
of the outer planet. 

The density in the system is adjusted such that in the 
relaxed initial state there is approximately 2.75 x 1O“^M0 
of material within the computational domain (see below). 

3.5. A few remarks on numerical issues 

The numerical method used is a staggered mesh, spatially 
second order finite difference method, where advection is 
based on the monotonic transport algorithm (van Leer 
1977). The code uses operator-splitting and is semi-second 
order in time. The computational details of the code that 
we employ, RH2D, have been described in general in Kley 
(1989), and specifically for planet calculations in Kley 
(1999). The use of a rotating coordinate system requires 
special treatment of the Coriolis terms to ensure angular 
momentum conservation (Kley 1998), an especially impor¬ 
tant point for the long-term calculations presented here. 

The viscous terms, including all necessary tensor com¬ 
ponents, are treated explicitly. To ensure stability in the 
gap region, where there are very strong gradients in the 
density, an artificial bulk viscosity has been added, with 
a coefficient Cart = 2. For a detailed discussion of the 
viscosity related issues and tests, see Kley (1999). 

The energy equation Eq. (1) is solved explicitly apply¬ 
ing operator-splitting. The heating and cooling term D—Q 
is treated as one sub-step in this procedure, see Gunther 
et al. (2004). The additional radiative diffusion part in the 
energy equation (5) is solved applying an implicit method 
to avoid possible time step limitations. To solve the result¬ 
ing matrix equations we use Successive Over Relaxation 
(SOR) with an adaptive relaxation parameter (Kley 1989). 

A larger mass ratio M/M, induces stronger torques 
and produces low densities in the gap region. To prevent 
numerical instabilitites caused by too large gradients we 
have found it preferable to work with a density floor, where 
the density cannot fall below a specified minimum value 
Smin- For our purpose we use a value of Emin = 10“® in 


dimensionless values, where the typical (initial) density is 
of 0(1). 



Fig. 3. The azimuthally averaged density profile for the relaxed 
configurations for three different masses of the outer planet: 
§2 = 1.0 X 10 “®,g2 = 3.5 X 10“® and 5.9 x 10“®. The planet 
is located at a 2 = 0.35 with a fixed semi-major axis, and is 
allowed to accrete. In these models the mass of the inner planet 
has been switched off. The density is given in dimensionless 
units. 


3.5.1. Planetary dynamics 

The motion of the planets is integrated using a fourth 
order Runge-Kutta integrator where the time step size is 
given by the hydrodynamical time step. While not the 
most accurate integrator for full N-body calculations, it 
is sufficient for our purposes. As a test we have run the 
pure 3-body problem of one star with two planets under 
exactly the same conditions as in the full hydrodynamical 
evolution and find that the relative total energy loss is less 
than 6 x 10“® in 1000 years (~ 2.5 x 10® time steps), which 
is equivalent to over 6,000 periods of the inner and over 
around 2,700 periods of the outer planet. 

The forces of the disk are taken into account in a first 
order approximation to reduce the computational effort. 
To avoid problems with under-resolved material close to 
the planet, a torque cutoff radius of rtorq is applied where 
material closer to the planet than rtorq does not contribute 
to the force acting on the planet. The problem of choos¬ 
ing the optimum value for rtorq is non-trivial. Using a 
cutoff radius prevents large unphysical variations of the 
forces acting on the planet. Very high resolution simula¬ 
tions (D’Angelo et al. 2002, 2004) show a nearly symmet¬ 
ric distribution of material close to the planet indicating 
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that those regions do not contribute too much to the total 
torque. As these regions cannot be resolved in our sim¬ 
ulations we have to use a torque cutoff instead. For all 
isothermal models we use rtorq = 0.5i?Hiii, where the Hill 
radius is given by 

i?Hin = a2 (1)'^". (7) 

Tests with an increased value of 0.75i?Hiii gave indistin¬ 
guishable results on the migration rate and eccentricity 
evolution. For the radiative models we use a value of 
J’torq = 0.8i?Hiii because the radiative cooling leads to 
more material around the planet. 

In calculating the gravitational potential of the disk 
and planet the vertical extent of the disk is taken into 
account, assuming a vertically isothermal structure, i.e. 
Gaussian density distribution. For the smoothing length 
of the potential we choose rpot = 0.877. For a l.OMjup 
planet this is equivalent to only about 0.37?Hin. 

4. Model computations 

Constructing our models for the formation of GJ 876 con¬ 
sists of two basic steps: 

I. Construct equilibrium models where the planets have 
fixed orbital elements, and let the disk evolve into a quasi¬ 
stationary state. 

II. Follow the subsequent evolution by ‘releasing’ the plan¬ 
ets, i.e. by taking into account the disk forces acting on 
the planet. 

This two-fold procedure is necessary to avoid evolutions 
which might be dominated by a transient adjustment of 
the disk, as the disk-planet equilibrium state is not known 
a priori. The disadvantage is that the equilibration phase 
may take many hundreds of orbits. Let us consider these 
two steps in turn. 

4.1. Relaxation towards equilibrium 

4.1.1. Switched off inner planet 

In the first part of step I the mass of the inner planet 
is switched off (reduced by a factor 10“®) and the outer 
planet has fixed orbital elements (02 = 0 . 35,62 = 0.01). 
We construct equilibrium models for different planet 
masses <72 because firstly, the observations do not yield 
definite masses of the planets due to the poorly con¬ 
strained inclination, and secondly we want to understand 
the planet-disk system in more general terms. If the sur¬ 
face density and all other variables are fixed at the outer 
boundary there corresponds (for a given viscosity and ver¬ 
tical height of the disk) one particular equilibrium state 
to each mass ratio. 

In Fig. 3 the azimuthally averaged E(r) profile is shown 
for the relaxed equilibrium disk states for three different 
masses q 2 of the outer planet starting from (72 = 1.0 x 
10“® (black solid line), over 3.5 x 10“® (gray dotted line) 
upto 5.9 X 10“® (dashed line), with the last value close 



Fig. 4. Gray scale plots of the surface density E for the re¬ 
laxed state with no inner planet for two different masses: top) 
(72 = 3.5 X 10“® and bottom) <72 = 5.9 x 10“®. Due to the 
higher planetary mass much stronger wave-like disturbances 
are created in the density, and the disk becomes eccentric with 
very small pattern speed in the inertial frame. 


to the minimum mass of the outer planet of the GJ 876 
system. One notices that for the two smaller planetary 
masses the outer disk is much more “quiescent” than for 
the large mass. Larger planetary masses apparently lead 
to a restructuring of the disk. The wave-like disturbances 
in the disk (seen as the oscillatory behavior in the curves) 
are significantly stronger for (72^5.2 x 10“®. 

The existence of the two different equilibrium states 
of the disk is further illustrated in Fig. 4 where we dis¬ 
play gray scale plots of the surface density S for the re- 
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Fig. 5. The time evolution of the semi-major axis and eccentricity for three models with different masses of the outer planet, 
and zero mass inner planet. 


laxed state with no inner planet for two different masses 
(g2 = 3.5 and 5.9 x 10“^). In both cases the small but 
non-vanishing eccentricity of the outer planet (e 2 = 0.01) 
leads to wave-like disturbances in the disk oscillating with 
the planetary period. But while for the lower mass case 
(92 = 3.5 X 10“^) the disk structure remains quite regu¬ 
lar, the second high mass case (92 = 5.9 x 10“^) shows 
a strongly disturbed disk which has gained a significant 
eccentricity of about 0.25 near the gap edge. The tran¬ 
sition from the non-eccentric state to the eccentric state 
which is here a function only of the planetary mass may 
depend also on the viscosity and temperature on the disk 
which we have held hxed. To check if this effect is in¬ 
duced by the non-zero eccentricity of the planet we have 
run a model with a vanishing eccentricity (e 2 = 0) and 
found the same behavior. We find that the eccentric disk 
is nearly stationary in the inertial frame. The mass flow 
onto the planet and through the gap is higher than the 
flow for the non-eccentric disk that prevailed for the lower 
mass planet. While it would be interesting, a further ex¬ 
ploration of these features is beyond the scope of this pa¬ 
per, and we leave it to a subsequent study. This disk ec¬ 
centricity may be closely related to the disk eccentricity 
induced by resonant interaction of the outer 3:1 eccentric 
Lindblad resonance of planet 2 with the disk, as discussed 
by Papaloizou et al. (2001), and to the problem of the in¬ 
ner disk in cataclysmic variables by Lubow (1991). This 
violent interaction of a massive planet with the disk ex¬ 
plains also some of the outlined numerical difficulties and 
requirements (such as density floor and artihcial viscos¬ 
ity)- 



Time [Yrs] 


Fig. 6. The time evolution of the periastron for three models 
with different masses of the outer planet. The thick dashed line 
is a fit corresponding to a 1 deg/yr shift. 

4.1.2. Evolution of the outer planet alone 

To explore the effect the disk alone has on the evolution of 
the orbital elements of the single (outer) planet, in partic¬ 
ular the influence on the precession rate, we present some 
calculations where the gravitational back-reaction of the 
disk is taken into account, while the inner planet is not 
present. For this purpose we chose the density unit such 
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that S = 1 (e.g. in Fig. 3) refers to a surface density of 
approximately 7500 g cm“^. This value refers to a total 
mass of 2.75Mjiip within the computational domain for 
the q 2 = 0.0059 case. 



Fig. 7. The relaxed azimuthally averaged density profile for 
two isothermal models with and without considering the inner 
planet (solid, dotted lines), and two radiative models includ¬ 
ing only heating and cooling (short-dashed), and additionally 
radiative diffusion (long-dashed). The mass of the inner planet 
is qi = 1.75 X 10~® and that of the outer 52 = 5.9 x 10“®. The 
gray short-dashed line is identical to the corresponding line in 
Fig. 3 (labeled 5.9). 

In Figs. 5 and 6 the time evolution of the orbital el¬ 
ements of the outer planet for three models with differ¬ 
ent planet masses are displayed. All models have been re¬ 
laxed before releasing the planet. For direct comparison 
the physical disk mass in all three cases has been scaled 
to be identical (« 2.75 x 10 “^Mq). With respect to the 
reference model (c3) with 52 = 0.0059, this implies for 
models (d3a) and (e4a) a density reduction factor of 0.98 
and 0.9, respectively (cp. Fig. 3). The migration rate of 
the planets depends primarily on the amount of mass near 
the 2:1 Lindblad resonance. We observe that the lowest 
{q 2 = 0.001) and highest {q 2 = 0.0059) mass planet have 
very similar initial migration rates, while the intermediate 
mass planet has a faster migration rate. The intermediate 
mass {q 2 = 3.5 x 10“^) has the largest rate because there is 
more mass very close to the location of the 2:1 resonance, 
which lies here at r = 0.56 (note the steep density gradient 
in Fig. 3). From the average density profile one might have 
suspected a smaller migration rate for the massive planet. 
However, due to the eccentric disk (Fig. 4), material gets 
close to the planet every orbit, which increases the mi¬ 
gration rate slightly. Thus, the non-monotonic variation 


of migration speed with planet mass is related to the dif¬ 
ferent equilibrium disk states, circular and eccentric. The 
right panel of Fig. 5 shows a temporary increase of the ec¬ 
centricity followed by a decline for the larger mass planet 
case, while the smaller masses both show declining eccen¬ 
tricities at nearly the same rates indicating the standard 
eccentricity damping in planet-disk interaction. This dif¬ 
ferent behavior of the eccentricity evolution is definitely 
caused by the eccentric disk state for the higher planet 
mass case. It remains to be studied if this effect may have 
some relation to the observed large eccentricities of the 
extrasolar planets. 

In Fig. 6 the evolution of the periastron is displayed for 
the three cases. For comparison the superimposed black 
dashed-line refers to a periastron advance of exactly l°/yr. 
This positive Wdisk ~ l°/yr (for the small and high mass 
planet) is solely due to the interaction of the disk with 
outer planet alone. The amount of planetary precession 
induced by the presence of the disk is determined primar¬ 
ily by the material being very close to the planet. Here, 
the behaviour of the low and high mass models are com¬ 
parable again, while for the intermediate mass, there is 
less material very close to the planet and the precession 
rate is smallest. 



.4 .6 .8 1 1.2 


r 


Fig. 8. The relaxed azimuthally averaged temperature profile 
for the isothermal (solid line) and two radiative models includ¬ 
ing only heating and cooling (short-dashed), and additionally 
radiative diffnsion (long-dashed). The mass of the inner planet 
is qi — 1.75 X 10“® and that of the outer 52 = 5.9 x 10~®. 


4.1.3. Switched on inner planet 

In the second step of the relaxation process the inner 
planet is included, using qi = 1.75 x 10“^. Additionally, we 
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use the full thermodynamics for some models. The effect 
that these additions have on the density and temperature 
distribution is displayed in Figs. 7 and 8. The surface den¬ 
sity is slightly increased upon including the inner planet 
because the additional torques tend to push the matter 
a bit further away from the central star. The presence of 
the outer planet is now seen clearly in the surface den¬ 
sity distribution because for these models the planet is 
not allowed to accrete material from its Roche-lobe. The 
subsequent accumulation of gas in the Roche-lobe appears 
as a spike near r = 0.35 in Fig. 7. The inclusion of radia¬ 
tive diffusion (Eq. 5) in the plane of the disk reduces the 
temperature in the region of the planet and leads to a 
larger mass accumulation in the Roche lobe of the planet 
than in the model which includes only heating and cool¬ 
ing (Eq. 1). The density distribution in the disk and gap 
region are not influenced too much by including radiative 
diffusion. 
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Fig. 9. The time evolution of the orbital elements (a, e, 9i, Arcj) 
for an isothermal model (h4) with qi = 1.75 x 10“® and q 2 = 
3.5 X 10“^ (test case with lower outer planet mass). 


Fig. 10. The time evolution of the orbital elements 
(a, e, 01 , Aro) for an isothermal model (h2) with gi = 1.75 x 
10“® and q 2 = 5.9 x 10“® (as GJ 876 edge on). 




Time [Yrs] 


Time [Yrs] 


In Fig. 8 the change in temperature due to the inclu¬ 
sion of radiative effects is displayed. Due to the relatively 
large density of disk and the high viscosity, the temper¬ 
ature in the radiative cases rise considerably above the 
isothermal case, and is given in the disk region by the equi¬ 
librium of heating and cooling D = Q. The disk thickness 
increases from H/r = 0.05 to about H/r = 0.15, for the 
full radiative models. Including radiative diffusion in the 
plane of the disk (long-dashed line) leads to higher temper¬ 
ature in the gap region. The full thermodynamic models 
are no longer eccentric even for the large planetary mass. 
The included radiation and larger H/r leads to a narrower 
gap and additional damping of the modes. We do not in¬ 


Fig. 11. The time evolution of the orbital elements 
(a, e, 01 , Aro) for an isothermal model (h9) with q\ = 2.10 x 
10~® and 52 = 7.1 X 10“® (as GJ 876 at 55deg inclination). 


vestigate at this point the detailed dependence of the disk 
thermodynamics on variations of the physical viscosity. 

These fully relaxed models including the inner planet 
serve as the starting point for the dynamical evolution of 
the planets. The periastron advance for the outer planet 
due to the inner planet having qi = 0.00175 and ai = 
0.20 alone, with no disk forces is found to be tbpianet ~ 
0.68°/yr. Thus, the disk and the planet generate a shift of 
similar magnitude. 
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Fig. 12. The time evolution of the orbital elements 
(a,e,8i, Aro) for a radiative model including heating and cool¬ 
ing (ml) with qi = 1.75 X 10“® and q 2 = 5.9 x 10“® (as GJ 876 
edge on). Note that at around t = 380 the outer planet leaves 
the computational grid and the results are not reliable any¬ 
more. 


a test case (model h4) with a lighter outer planet (gi = 
0.00175,(72 = 0.0035), a case (h2) resembling the edge- 
on {qi = 0.00175, (72 = 0.0059) system, and finally the 
i = 55° {qi = 0.0021,(72 = 0.0071) case (h9), respectively. 
In all cases the outer planet captures the inner one into 
a 2:1 resonance. Although the orbits are initially quite 
far from 2:1 mean-motion commensurability, the resonant 
angle 9i is already in libration, with amplitude as large as 
~ 90°. However, this libration of 9i has little dynamical 
consequence and the increase in ei is slow, until the mean 
motions approach the 2:1 commensurability. There is a 
delay in the capture of Aw (and 02) into resonance, and 
the libration amplitude of 0i is smaller than that of Aw. 
Note that in the first test-case for smaller planet mass, 
when the initial disk state was not eccentric, the capture 
of Aw occurs faster and the libration amplitudes of both 
9 1 and Aw are much smaller. 

The fact that 9i (and even 02 and Aw) can be librat- 
ing when the orbital mean-motions are far from the 2:1 
commensurability is due to the l/cj dependence of the 
resonance-induced retrograde precession of Wj at small 
ej. The relative timing of the capture of 0i and 02 into 
resonance is affected by the masses and initial eccentric¬ 
ities of the planets. We shall discuss these points further 
in Sect. 5. 


Table 2. Parameters of the models used to study the evolution 
of a planetary system under the action of an outer disk. Given 
are a model name, the used thermodynamics (TD), isothermal 
(iso) or radiative (rad), the type of outer boundary condition 
(obc) (open or closed), the mass-ratio of the inner planet (qi), 
the mass ratio or evolution of the outer planet ( 92 ), the fig¬ 
ure in which the model is displayed. The displayed radiative 
models include heating and cooling only as the inclusion of the 
diffusion does not change results significantly 


Name 

TD 

obc 

91 

[10-3] 

0 ^ 

1 

to 

Fig. 

h4 

iso 

closed 

1.75 

3.50 

9 

h2 

iso 

closed 

1.75 

5.9 

10 

h9 

iso 

closed 

2.10 

7.10 

11 

ml 

rad 

closed 

1.75 

5.9 

12 

k3b 

iso 

Open 

1.75 

5.9 

13 

h8a 

iso 

open 

2.10 

5.9^ 6.9 

14 

m2 

rad 

open 

1.75 

5.9 

15 


4.2. Evolving planets 

Having constructed several equilibrium models holding 
the orbital elements of the planets fixed, we now release 
the planets and follow the evolution of their orbital ele¬ 
ments. For an overview Table 2 lists the models and their 
most important parameters. 

In the hrst instance we keep the damped and reflecting 
boundary conditions at the outer boundary, i.e. we model 
the situation where the disk remains in full contact with 
the planet during the evolution. In Figs. 9 to 11 we display 
the evolution of the semi-major axes, the eccentricities, 
and the two resonant angles 0i and Aw for three models. 


Once both 0i and Aw are librating about 0°, the ec¬ 
centricities rise rapidly. In cases h2 and h9, applicable to 
GJ 876, the eccentricity of the inner planet rises above 
0.4, which clearly exceeds the upper limit of ^ 0.31 for 
GJ 876 (see Fig. 1). 

The evolution with the inclusion of heating and cool¬ 
ing (ml) is displayed in Fig. 12. Again the system is 
caught in a 2:1 resonance, and the eccentricities rise to 
high values. Here, due to the smoother initial state, the 
alignment of the resonant angles is much stronger and the 
libration amplitudes for both angles are reduced strongly. 
The evolution beyond t = 380 yr is unreliable because 
both planets are outside the computational domain by 
then. The models with radiative diffusion do not yield 
any different results from models with only heating and 
cooling. The equilibrium conhgurations as shown in Figs. 7 
and 8 above are very similar already. 

In summary, there exist two different types of evolu¬ 
tions of the orbital elements depending on the state of 
disk. For a nearly circular disk (models h4, ml) the de¬ 
lay in the capture of Aw (or equivalently 02 ) allows 62 
to be initially damped. After capture in this second reso¬ 
nance 62 rises rapidly and increases far beyond the upper 
limit inferred from the observations. The libration of the 
resonant angles after capture is very small. This result 
confirms earlier findings of hydrodynamical evolutions of 
resonant planets (Kley et al. 2004). For an eccentric disk 
state, i.e. for larger outer planet masses, capture of 02 is 
more delayed and the libration of the two resonant an¬ 
gles remains quite large. Nevertheless, the eccentricities 
quickly exceed the upper limits for GJ 876. 
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Fig. 13. The time evolution of the orbital elements 
(a, e,8i, A-cj) for an isothermal model (k3b) with qi = 1.75 x 
10“^ and 52 = 5.9 x 10“®, including disk dispersal. 


Fig. 15. The time evolution of the orbital elements 
(a, e, 01, Azu) for a radiative model (m2) with qi = 1.75 x 10“® 
and 52 = 5.9 x 10“®, including disk dispersal. 




Fig. 14. The time evolution of the orbital elements 
(a,e,di, Azj) for an isothermal model (h 8 a) with 51 = 1.75 x 
10~® and a growing outer planet 52 = 5.9 ^ 6.9 x 10“^. The 
model includes disk dispersal. 


beyond i?max) to 1/10 of the value of the adjacent inner 
gridcell. This results in a mass loss rate which is some¬ 
what higher than physically plausible, but allows us to 
halt the migration in a computationally feasible time. For 
the isothermal models (k3b) and (h 8 ) we find a mass- 
half-emptying time of about 300 yrs and for the radiative 
model 100 yrs. 

In Figs. 13 and 14 two isothermal cases are presented; 
in the first (k3b), the mass of the outer planet remains 
constant, while in the second model (h 8 a) the planet is 
allowed to accrete from its surroundings and to grow in 
mass. In the second case the mass of the inner planet was 
chosen higher 51 = 0.0021 to compensate for the higher 
final mass of the outer planet. Again, capture into the 2:1 
resonance occurs, but the eccentricities do not rise to such 
high values. The loss of mass from the disk turns the mi¬ 
gration off before the eccentricities get large. The align¬ 
ment of the resonant angles occurs on longer timescales 
than before. Similar behavior can be seen in the evolution 
of the radiative model (m2) displayed in Fig. 15. Again the 
alignment of the resonant angles occurs faster and libra- 
tion of the resonant angles is smaller than in the isother¬ 
mal models due to the non-eccentric disk state. 


4.2.1. Modeling disk dispersal 

The next set of models starts with the same conditions 
as the previous ones but allows for material to leave the 
outer boundary; this could occur if the outer disk is pho- 
toevaporated, for example. Depending on the maximum 
outflow velocity allowed, the mass in the disk can be re¬ 
duced slowly or very efficiently. For our purposes we set 
the velocity in the outer radial ghostcell (the gridcell just 


5. Interpretations 

In this section we use three-body integrations and analytic 
theory to interpret the hydrodynamic results presented in 
Sect. 4. The three-body integrations were performed using 
the symplectic integrator SyMBA (Duncan et al. 1998), 
modified to include forced migration and apsidal preces¬ 
sion and to have input and output in Jacobi coordinates 
(Lee & Peale 2002; Lee 2004). Although three-body inte¬ 
grations do not model the interactions between the disk 
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and the planets self-consistently, they allow different ef¬ 
fects to be included separately to assess their importance. 
In particular, we show that the disk-induced apsidal pre¬ 
cession, which has not been included in previous three- 
body integrations, does not cause an additional delay in 
the capture of 02 and Anj, and that the additional delay 
in the capture of A.w in, e.g., the isothermal models h2 
(Fig. 10) and h9 (Fig. 11) is due to the relatively large 
initial eccentricities induced by the eccentric disk. 






Fig. 16. Time evolution of the orbital elements (a, e, Ara) 
for the baseline three-body model with qi = 1.75 x 10“® and 
q 2 = 5.9 X 10~®. The planets are initially on coplanar circular 
orbits. The outer planet is forced to migrate inward at the 
rate 02/02 = —4.8 x 10“^ (0.35 AU/ 02 )®^^ yr“^, and there is 
no eccentricity damping or additional apsidal precession. 


Figure 16 shows the evolution of the semimajor axes, 
eccentricities, and resonant angles 9i and Atu for the 
baseline three-body model with qi = 1.75 x 10“^ and 
q 2 = 5.9 X 10“^. The outer planet is forced to migrate in¬ 
ward at the rate 02/02 = —4.8x 10“^ (0.35 AU/ 02 )^/^ yr“^ 
(which matches the migration rate at 02 = 0.35 AU for the 
h2 model in Fig. 10), and there is no eccentricity damp¬ 
ing or additional apsidal precession. The planets are ini¬ 
tially on coplanar circular orbits, with oi = 0.2 AU and 
02 = 0.4 AU. The circular orbits and larger 02 compared 
to the hydrodynamical simulations (where 02 = 0.35 AU) 
reduce the initial eccentricity variations. As in the hy¬ 
drodynamical simulations, the resonant angle 61 is al¬ 
ready librating about 0° at the start of the evolution, even 
though the planets are started even further from the 2:1 
mean-motion commensurability than in the hydrodynam¬ 
ical simulations. With the small initial eccentricity varia¬ 
tions and the use of Jacobi coordinates, it can be seen in 
Fig. 16 that Atu (and 62 ) are first captured into libration 


about 180° at t « 270 yr and that Avu spends more time 
around 180° even before the capture. The system passes 
smoothly over to the configuration with both 0i and Azu 
librating about 0° when ei ^ 0.1 (see also Lee 2004). The 
small libration amplitudes of the configuration with both 
61 and Azu librating about 0° are similar to those found 
in the radiative model ml (Fig. 12). 

For two planets orbiting a star in coplanar orbits, the 
equations of motion for the periapse longitude and eccen¬ 
tricity are 


d'uj j 

dt 


dej 

dt 




9$ 


MjCj^GA4:^cLj dtcj 


T dUsec^j T ff^diskj' 




MjCj ^GM^aj dzuj 


(1 - e]) - 


1 


Mj e j ^ Glvl^cij d Xj 


Cdiskjj 5 


( 8 ) 


(9) 


where the disturbing potential 


GM1M2 , 

4 > —- (/?(/ 3 , ei, 62, 6*1, 62), 

02 


( 10 ) 


the ratio 01/02 = /3, and is a function of the indi¬ 
cated variables if we consider only the 2:1 resonant terms 
and neglect terms of order [(Mi + M 2 )/M*]^ and higher 
(e.g. Yoder & Peale 1981; Lee & Peale 2002; Beauge & 
Michtchenko 2003). The term Wdiskj represents the posi¬ 
tive apsidal precession induced by the disk, zusecj the di¬ 
rect secular effect from the other planet, and edisk.j the ec¬ 
centricity damping induced by the disk. The disk-induced 
variations are not included in Fig. 16, and we neglect them 
in our discussion for the moment. 

When the eccentricities are very small, the relevant 
terms in ip are C'i(/J)ei cos^i and C 2 (/J)e 2 cos 02 , and we 
have 


dzui/dt = (3q2niCiei^ cos9i + vusec,i, (11) 

dw2/dt = gin2C'2e^^ COS02 + ff7sec,2, (12) 

and 

dei/dt = — f 3 q 2 niGi sin 9i, (13) 

de2/dt = —(jfin2C'2 sin02, (14) 


where rij are the mean motions and Ci{(3) « —1.19 and 
C'2(/3) ~ -1-0.43 for [3 « 2“^/^ (Yoder & Peale 1981). Stable 
simultaneous librations of both 9i and 02 (or equivalently 
Azu) require that the longitudes of the periapses regress 
at the same rate on average and that the eccentricities do 
not change in the absence of continued migration. Since 
Cl < 0 and C 2 > 0, the only way these requirements are 
satisfied is for 9i to librate about 0° and 02 and Azu to 
librate about 180°. This anti-aligned configuration is what 
we observe with the lo-Europa pair of Jupiter’s satellites 
and is shown in Fig. 16 just after the capture of Azu. 

The relative timing of the capture of 0i and 02 into 
resonance needs some clarification. Due to the 1/ei de¬ 
pendence of the resonance-induced precession of vu\ at 
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small ei (Eq. (11)), a small initial value of ei leads to an 
extremely mobile longitude of periapse wi, such that the 
large mass of the outer planet can induce sufhcient ret¬ 
rograde motion of w\ to cause 6 *i to already be librating 
at the start of the evolution when the orbital mean mo¬ 
tions are far from the 2:1 commensurability. This initial 
libration is evident in the three-body model shown in Fig. 
16 and in the hydrodynamical simulations shown in Sect. 
4.2. It should be noted, however, that this libration of B\ 
has little dynamical consequence and the increase in ei 
(see next paragraph) is slow, until the mean motions ap¬ 
proach the 2:1 commensurability. In spite of the relatively 
small initial value of ei in Fig. 16 and the same 1/e de¬ 
pendence of the resonance-induced precession (Eq. (12)), 
02 and hence Atu are not librating initially, although they 
do spend more time near 180°. This is due to the inner 
planet’s lower mass being insufficient to perturb the more 
massive outer planet’s vui enough to cause 6*2 to librate 
that far from the 2:1 commensurability of the mean mo¬ 
tions. 

The reason the eccentricities increase with continued 
migration can be understood as follows. Continued inward 
migration of the outer planet while the system is within 
the resonances means ai is consistently slightly smaller 
that it would have been without the migration. The in¬ 
creased value of ni therefore means A 2 is slightly larger at 
any instant than it would be without the migration. The 
argument of the sine in Eq. (13), di = 2 X 2 — Ai — tui, is 
then slightly greater than 0° on average. Since Ci < 0, 
dei/dt > 0 and the eccentricity must grow. A similar ar¬ 
gument applies to dei/dt if 62 is also librating about 180° 
(see Eq. (14)). It can be shown from the equations for the 
evolution of the orbital energy and angular momentum 
that at least one, but not necessarily both, of the eccen¬ 
tricities must be increasing at any given time for continued 
migration within the resonance, if < (3 -I- ef)/4 for the 
2:1 resonance and there is no eccentricity damping (see 
Eq. (AlO) of Murray et al. 2002). 

For Ml/Ml « 0.31, the resonant interaction is no 
longer dominated by the lowest order resonant terms when 
Cl ^ 0.04, and ei changes from increasing to decreasing 
with continued migration until 62 ~ 0 when ei « 0.1. This 
change in the evolution of ei is not obvious in Fig. 16 
because ei ^ 0.002 during this phase; see Fig. 5 of Lee 
(2004). Then the libration center of 9i and Atu changes 
from 180° to 0 ° (with a slight phase shift due to the contin¬ 
ued migration), and both eccentricities continue to grow 
because of the slight phase shift in the arguments of the 
terms involving linear combinations of 9i and Bi. However, 
because there are now many terms in the disturbing po¬ 
tential that contribute to the resonant interaction, a sim¬ 
ple demonstration of the co-precession of the periapses 
with both 6*1 and 9i librating about 0 ° and of the in¬ 
crease in both eccentricities with continued migration is 
elusive. Continued migration causes the eccentricities to 
quickly exceed those observed without a large eccentricity 
damping as discussed above. Although the hydrodynam¬ 
ical simulations self-consistently produce Cdisk < 0 , the 


magnitude of this damping is far less than that necessary 
to maintain the small observed eccentricities in the GJ 876 
system during continued migration. 



Fig. 17. Time evolution of the orbital eccentricities for the 
three-body models with additional apsidal precession d7disk,2 = 
p (0.35 AU/a 2 )®^'^ ° yr“^ and p = 1 and 10, compared to the 
case without additional apsidal precession (p = 0) from Fig. 
16. 


5.1. Disk-Induced Apsidal Precession 

It is shown in Fig. 6 that the disk induces a prograde ap¬ 
sidal precession of the outer planet of about l°yr~^ for 
qi = 5.9 X 10“^. This prograde apsidal precession is pri¬ 
marily due to the axisymmetric component of the disk 
potential. For a planet with orbital semimajor axis a and 
mean motion n, an outer disk with surface mass density 
S(r) = So(ro/r)* at ri < r < r 2 (where ri > a) induces 
a prograde apsidal precession (Ward 1981) 


tbdisk = 27rn (^) [Wk(ri/a) - Wk{ri/a )], (15) 

where 


WMU) = y: 


i=l 


m 


22<?(H)2J Vr 


(- 


2^+fc-l 
) • ( 16 ) 


If we approximate the azimuthally averaged density profile 
for qi = 5.9 x 10“^ in Fig. 3 as a flat profile (fc = 0) with 
So = 7500gcm“2 between ri = 0.6 AU and ri = 1.2 AU, 
we find 'd 7 disk ,2 ~ 0?7yr“^, which is in reasonable agree¬ 
ment with the measured value. The disk-induced preces¬ 
sion rate can be increased by increasing the disk mass or 
decreasing the gap width (ri — a)/a, but the gap width 
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is determined mainly by the mass ratio q and is not very 
sensitive to either the viscosity parameter a or the disk 
thickness H/r (e.g. Kley 1999; Bryden et al. 1999; Varniere 
et al. 2004). Increasing a and/or H/r primarily makes the 
gap shallower and not so deep. 

We have performed three-body integrations similar to 
that in Fig. 16, but with additional apsidal precession 
^i 7 disk ,2 = p (0.35 AU/a 2 )^/^ and p = 1 and 10. The 
much smaller disk-induced precession of the inner planet’s 
orbit was neglected. Except for a larger libration ampli¬ 
tude for Atz 7 in the p = 10 case, the evolution of 9i and 
Atu is similar to that shown in Fig. 16 for p = 0, and 
there is no additional delay in the capture of Atu. Figure 
17 shows the evolution of the eccentricities for p = 0, 1, 
and 10. At a given time, the additional positive apsidal 
precession results in an increase in ei and a decrease in 
62 , but the effect is small for ei < 0.31 (the upper limit 
for GJ 876; Fig. 1) even for p = 10. The additional apsidal 
precession and forced migration rates used in the three- 
body integrations are proportional to 02 and hence the 
mean motion n 2 - A decrease in the positive, disk-induced 
precession rate due to, e.g., disk dispersal would cause the 
eccentricities to adjust to values appropriate for the pre¬ 
cession rate at a given time and approach the p = 0 case 
without disk-induced precession. 

The analytic theory developed by Yoder & Peale 
(1981) for the 2:1 resonances of lo and Europa and the 
Laplace resonance takes into account the apsidal preces¬ 
sion induced by the oblateness of Jupiter (which is largest 
for the innermost satellite lo), and it can be adapted to 
understand the effects of the disk-induced apsidal pre¬ 
cession (which is larger for the outer planet). As the or¬ 
bits converge toward the 2:1 mean-motion commensura- 
bility (i.e., as 2 n 2 — ni < 0 increases), the resonance 
6*1 = 2 A 2 — Al — 1371 would be encountered before the 
resonance 6*2 = 2 A 2 — Ai — tz 72 if the apsidal precession 
is dominated by Wdiskj (since 9j « 2 n 2 — ni — zuj and 
d 7 disk ,2 i77disk,i > 0). However, because the resonance- 
induced retrograde apsidal precession is proportional to 
1 /ej and much larger in magnitude than the disk-induced 
prograde precession for small e^, 9\ and 9^ (and hence 
Atu) can be captured into libration in a sequence that 
differs little in order or timing from the case where there 
is no disk-induced precession. When we include Wdiskj' in 
Eqs. (11) and (12) for small e^, stable retrograde preces¬ 
sions of both orbits still require ( 0 i) = 0 ° and (^ 2 ) = 180°, 
and the requirement that the orbits process at the same 
rate on average ((tiji) = ( 1772 )) implies the following rela¬ 
tionship between the forced eccentricities: 


62/61 = qin^C^/X-^qinxCx J- ei(d7sec.2 - djsec.i) 

-t -61 ('d7disk,2 ^disk,l)]- (^^) 

Since 'ri 7 disk ,2 — ^diskp > 0, the disk-induced apsidal pre¬ 
cession reduces 62 / 61 . The decrease in 62 and increase in 
61 seen in Fig. 17 are consistent with this trend, but it 
should be noted that the 2:1 resonance configurations with 
M 1 /M 2 « 0.3 and Ci ^ 0.04 are no longer dominated by 


the lowest order resonant terms and those with Ci ^ 0.1 
also have 6*2 and Atu librating about 0 ° instead of 180° 
(see above). Although 62/61 is fixed by the resonance con¬ 
ditions, a simple expression like Eq. (17) does not exist for 
the larger eccentricities. 
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Fig. 18. Same as Fig. 16, but for the three-body model with 
the following initial conditions: ei = 0.01, 62 = 0.05, and the 
orbits are antialigned, with the inner planet at periapse and 
the outer planet at apoapse. 


5.2. Initial Eccentricities 

The isothermal model h2 (Fig. 10) shows large initial ec¬ 
centricity variations induced by the eccentric disk, with 
both 61 and 62 around 0.05. The large eccentricities (in 
particular 62 ) mean that the planets cannot be captured 
into the 2:1 resonance configuration with 9i librating 
about 0° and 92 librating about 180° (see Fig. 16). In 
Fig. 18 we show the results of a three-body integration 
similar to the baseline model in Fig. 16, but with the fol¬ 
lowing initial conditions: 61 = 0.01, 62 = 0.05, and the 
orbits are antialigned, with the inner planet at periapse 
and the outer planet at apoapse. These initial conditions 
were chosen so that the eccentricity variations are similar 
to those in the model h2 when 02 « 0.35 AU. The larger 
value of 61 has prevented the initial libration of 9i that is 
seen in Fig. 16. But 9i is captured into libration about 0° 
at about the time {t « 300 yr) when 02 « 0.35 AU, and its 
libration amplitude is similar to that in model h2. There 
is also a delay in the capture of Aw (and 02 ) into libration 
to t fn 600 yr and, as expected, Aw is captured directly 
into libration about 0°. The libration amplitude of Aw is 
smaller than that in model h 2 , possibly because of contin¬ 
uing interaction with an eccentric disk in the latter case. 
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The relative timing of the captures can have no effect on 
the necessity of large eccentricity damping if there is ex¬ 
tensive migration within the resonances, since the ratio of 
the eccentricities is determined uniquely by the resonance 
conditions once all the resonance variables are librating. 

6. Conclusion 

We have modeled the evolution of the GJ 876 system by 
performing two-dimensional hydrodynamical simulations 
of a disk with embedded planets. We conhrm previous 
work showing that interactions between the outer planet 
and a gas disk can drive two initially non-resonant planets 
into resonance, and in particular we have shown that the 
induced migration of the GJ 876 planets results in the cap¬ 
ture of the 2:1 resonance variables, 9i, 02 and Aw into the 
observed libration about 0°. The precession of the outer 
planet’s periapse longitude induced by the disk was shown 
not to increase the delay of capture of 62 into resonance, 
but a high value of 62 induced by an eccentric disk did 
increase the delay. Later capture of 62 has no effect on the 
eccentricity damping problem as the ratio of the forced 
eccentricities is fixed once both 9i and 6*2 are librating. 

The “isothermal” simulations with fixed disk tempera¬ 
ture structure have shown that more massive planets with 
a planet-to-star mass ratio oi q > 5.2 x 10“^ are able to 
perturb the disk sufficiently to make it eccentric (even if 
the planetary orbit is circular). This effect may be caused 
by an eccentric instability driven by an interaction of the 
m = 2 mode at the outer eccentric 3:1 Lindblad resonance 
with a slightly eccentric disk (see Papaloizou et al. 2001). 
For smaller masses in the isothermal models and for all 
masses in the radiative models, the planet-disk interac¬ 
tion produces the typical gap and spiral arms, but the 
disk remains otherwise circular. 

The simulations also conhrm that for non-eccentric, 
circular disks the eccentricity growth of the outer planet 
is suppressed until 02 is locked into libration about 0 ° (see 
models h4, ml, m2). However, 02 and Aw are always even¬ 
tually captured into libration, and the subsequent increase 
in the eccentricities of both planets past the observational 
upper limits in GJ 876 (ei « 0.31, 62 « 0.05) occurs on 
a time scale shorter than the expected viscous time scale 
of the disk. The final eccentricities of both planets are 
then substantially larger than those seen in GJ 876, a re¬ 
sult which was found already in previous hydrodynamic 
(Papaloizou 2003; Kley et al. 2004) and three-body simu¬ 
lations (e.g. Lee & Peale 2002). 

In addition to causing an increased delay in the capture 
of 02 and Aw, eccentric disks (models h 2 , h9) lead to 
significantly larger libration amplitudes for Aw (« 50°) 
than for the circular disk case. Interestingly, it appears 
that in the case of GJ 876 the libration of Aw is in fact 
larger (« 34° for i = 90°) than that of 0i (« 7° for i = 
90°), as shown recently by Laughlin et al. (2004). 

Lee & Peale (2002) have shown that if sufficient ec¬ 
centricity damping ( 02/62 « 10002 / 02 ) is applied during 
migration, the observed configuration of GJ 876 can be 


maintained for arbitrary migration times. We have up¬ 
dated the amount of eccentricity damping required using 
the updated dynamical fits by Laughlin et al. (2004). For 
the best fits with coplanar inclination i ^ 35°, ei ^ 0.31, 
62 ^ 0.05, and 62/62 must exceed « 4002 / 02 . However, 
no such rapid eccentricity damping mechanism is presently 
known. The hydrodynamic simulations typically give com¬ 
parable time scales for eccentricity damping and semi¬ 
major axis decrease. Some current theories of planet-disk 
interactions indicate that eccentricity driving should oc¬ 
cur (Ogilvie & Lubow 2003; Goldreich & Sari 2003). Such 
eccentricity behavior would not be consistent with the ob¬ 
served state of GJ 876 if there were extensive migration 
after capture into the resonances. 

If the disk is removed soon after the planets are cap¬ 
tured into resonance, as in Figs. 13 to 15, the driving dis¬ 
appears and the final eccentricities are smaller, like those 
observed. However, the disk must vanish before the post¬ 
capture orbits shrink by ~ 10%, similar to the result found 
in 3-body integrations by Lee & Peale (2002). We have 
used in the present simulations a high viscosity and high 
surface density to increase the migration rates, and to be 
able to perform the simulations in a reasonable amount of 
computer time. In turn this implies a need for rapid disk 
dispersal to halt migration. If more realistic values were 
used for the viscosity and the disk mass, the necessary disk 
dispersal time scale could be much more extended. What 
appears to be somewhat fine tuned in the present simu¬ 
lations would be perhaps more reasonable. Our scenario 
only requires that capture occurred in the hnal formation 
phase of GJ 876, and that the planets did not migrate over 
a very long distance while locked in the resonances. 

Additional effects that might influence the eccentricity 
damping in these type of simulations are possibly three- 
dimensional. Fully 3D-simulations of circular planets on a 
fixed orbit (D’Angelo et al. 2003b; Bate et al. 2003) have 
shown that the spiral structures are not so clear and ma¬ 
terial may enter the Roche lobe from above the midplane, 
an effect which may alter the eccentricity damping prop¬ 
erties. The inclusion of magnetic helds might lead to a 
magnetic coupling of the planetary field with that of the 
disk, which will have an influence on the eccentricity evo¬ 
lution. But if the consideration of these and perhaps other 
processes fail to produce sufficient eccentricity damping, 
the elimination of the disk shortly after the capture of the 
GJ 876 planets into the 2:1 resonances is a possible, albeit 
perhaps unsatisfying means of accounting for the observed 
low-eccentricity configuration. 
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